# ============================================================
# Project Title: Running the hedonic treadmill: Lateral 
# insecurity and foreign real estate investment
# Replication Script: Observational Results
# Date: 02/07/2026
# ============================================================

# ============================================================
# 0. Setup: clear environment, load packages
# ============================================================
rm(list = ls())

# Package loader (install if missing)
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman")

pacman::p_load(
  dplyr, haven, janitor, sjlabelled, jtools, interplot, ggplot2, gridExtra, grid,
  network, igraph, glmmTMB, intergraph, tidyverse, readr, mediation, DiagrammeR, cobalt,
  psych, webshot2, gt, foreign, car, showtext, ggeffects, extrafont, tibble,
  tidycensus, stringr, forcats, Hmisc, tidyr, patchwork, ggpubr, broom.mixed
)

# Optional: make output easier to read
options(stringsAsFactors = FALSE)

# ============================================================
# 1. Paths + load CCES data
# ============================================================

# ---- Path handling ----
# For replication, prefer a PROJECT ROOT with relative paths.
# If you keep setwd(), at least fail loudly if path doesn't exist.
project_dir <- "~/Desktop/Replication/"
if (!dir.exists(project_dir)) stop("Project directory not found: ", project_dir)
setwd(project_dir)

# ---- Load data ----
cces_data <- read_dta("unmatched.dta") %>%
  haven::as_factor() %>%
  mutate(rowid = row_number())

# Optional: export to CSV for transparency / cross-software use
write.csv(cces_data, "converted_from_dta.csv", row.names = FALSE)

# ============================================================
# 2. Clean + construct variables (individual-level)
# ============================================================

cces_data <- cces_data %>%
  mutate(
    # ZIP: 5-digit padded
    inputzip = stringr::str_pad(as.character(inputzip), width = 5, side = "left", pad = "0"),
    
    # Urbanicity collapse
    urbancity = forcats::fct_collapse(
      urbancity,
      "Urban" = "City",
      "Suburb" = "Suburb",
      "Town" = "Town",
      "Rural" = "Rural area",
      "Other" = c("Other", "skipped", "not asked")
    ),
    
    # Immigration attitudes index components
    imm_att1 = ifelse(CC24_323a == "Support", 1, 0),
    imm_att2 = ifelse(CC24_323b == "Oppose", 1, 0),
    imm_att3 = ifelse(CC24_323c == "Oppose", 1, 0),
    imm_att4 = ifelse(CC24_323d == "Support", 1, 0),
    
    # Index (NOTE: you currently use att1, att3, att4 only)
    imm_att_index = rowMeans(dplyr::across(c(imm_att1, imm_att3, imm_att4)), na.rm = TRUE),
    
    # Economic perceptions
    econ_pers = dplyr::case_when(
      CC24_301 == "Gotten much better" ~  2,
      CC24_301 == "Gotten somewhat better" ~  1,
      CC24_301 == "Stayed about the same" ~  0,
      CC24_301 == "Gotten somewhat worse" ~ -1,
      CC24_301 == "Gotten much worse" ~ -2,
      TRUE ~ NA_real_
    ),
    
    # Belief that immigration causes inflation
    imm_inflation = dplyr::case_when(
      UCR415_9 == "selected" ~ 1,
      UCR415_9 == "not selected" ~ 0,
      TRUE ~ NA_real_
    ),
    
    foreign_inflation = dplyr::case_when(
      UCR415_6 == "selected" ~ 1,
      UCR415_6 == "not selected" ~ 0,
      TRUE ~ NA_real_
    ),
    
    # Perceived household income change
    income_change = dplyr::case_when(
      CC24_302 == "Increased a lot" ~  2,
      CC24_302 == "Increased somewhat" ~  1,
      CC24_302 == "Stayed about the same" ~  0,
      CC24_302 == "Decreased somewhat" ~ -1,
      CC24_302 == "Decreased a lot" ~ -2,
      TRUE ~ NA_real_
    ),
    
    # Perceived price change
    price_change = dplyr::case_when(
      CC24_303 == "Decreased a lot" ~ -2,
      CC24_303 == "Decreased somewhat" ~ -1,
      CC24_303 == "Stayed about the same" ~  0,
      CC24_303 == "Increased somewhat" ~  1,
      CC24_303 == "Increased a lot" ~  2,
      TRUE ~ NA_real_
    ),
    
    # Employment: employed vs not
    employed = dplyr::case_when(
      employ %in% c("Full-time", "Part-time") ~ 1,
      employ %in% c("Temporarily laid off", "Unemployed", "Retired",
                    "Permanently disabled", "Homemaker", "Student") ~ 0,
      TRUE ~ NA_real_
    ),
    
    # Racial resentment
    rr1 = dplyr::recode(CC24_441a,
                        "Strongly disagree" = 1, "Somewhat disagree" = 2,
                        "Neither agree nor disagree" = 3, "Somewhat agree" = 4, "Strongly agree" = 5),
    rr2 = 6 - dplyr::recode(CC24_441b,
                            "Strongly disagree" = 1, "Somewhat disagree" = 2,
                            "Neither agree nor disagree" = 3, "Somewhat agree" = 4, "Strongly agree" = 5),
    racial_resentment = rowMeans(dplyr::across(c(rr1, rr2)), na.rm = TRUE),
    
    # Age and bins
    age = 2024 - birthyr,
    age_group = cut(age, breaks = c(18, 29, 39, 49, 59, Inf), labels = 1:5, right = TRUE),
    
    # Education (ordinal numeric)
    educ_cont = dplyr::recode(educ,
                              "No HS" = 1, "High school graduate" = 2, "Some college" = 3,
                              "2-year" = 4, "4-year" = 5, "Post-grad" = 6),
    
    # Gender binary (woman = 1)
    gender_bin = ifelse(gender4 == "Woman", 1, 0),
    
    # Income midpoint coding
    income_cont = dplyr::case_when(
      faminc_new == "Less than $10,000" ~ 5000,
      faminc_new == "$10,000 - $19,999" ~ 15000,
      faminc_new == "$20,000 - $29,999" ~ 25000,
      faminc_new == "$30,000 - $39,999" ~ 35000,
      faminc_new == "$40,000 - $49,999" ~ 45000,
      faminc_new == "$50,000 - $59,999" ~ 55000,
      faminc_new == "$60,000 - $69,999" ~ 65000,
      faminc_new == "$70,000 - $79,999" ~ 75000,
      faminc_new == "$80,000 - $99,999" ~ 90000,
      faminc_new == "$100,000 - $119,999" ~ 110000,
      faminc_new == "$120,000 - $149,999" ~ 135000,
      faminc_new == "$150,000 - $199,999" ~ 175000,
      faminc_new == "$200,000 - $249,999" ~ 225000,
      faminc_new == "$250,000 - $349,999" ~ 300000,
      faminc_new == "$350,000 - $499,999" ~ 425000,
      faminc_new == "$500,000 or more" ~ 500000,
      TRUE ~ NA_real_
    ),
    
    # Party ID
    party_id = dplyr::recode(pid7,
                             "Strong Republican" = 1, "Not very strong Republican" = 2,
                             "Lean Republican" = 3, "Independent" = 4,
                             "Lean Democrat" = 5, "Not very strong Democrat" = 6,
                             "Strong Democrat" = 7),
    
    # Demographics
    homeowner = ifelse(ownhome == "Own", 1, 0),
    white = ifelse(race == "White", 1, 0),
    race_black = ifelse(race == "Black", 1, 0),
    race_asian = ifelse(race == "Asian", 1, 0),
    race_hispanic = ifelse(race == "Hispanic", 1, 0),
    race_other = ifelse(!(race %in% c("White", "Black", "Asian", "Hispanic")), 1, 0),
    
    # Lateral insecurity + FREI attitudes
    lat_ins = dplyr::recode(UCR420,
                            "Not difficult at all" = 1, "A little difficult" = 2,
                            "Fairly difficult" = 3, "Difficult" = 4, "Extremely difficult" = 5),
    
    frei_att1 = dplyr::recode(UCR421,
                              "Strongly agree" = 3, "Agree" = 2, "Slightly agree" = 1,
                              "Neither agree nor disagree" = 0, "Slightly disagree" = -1,
                              "Disagree" = -2, "Strongly disagree" = -3),
    frei_att2 = dplyr::recode(UCR422,
                              "Strongly agree" = 3, "Agree" = 2, "Slightly agree" = 1,
                              "Neither agree nor disagree" = 0, "Slightly disagree" = -1,
                              "Disagree" = -2, "Strongly disagree" = -3),
    frei_att3 = dplyr::recode(UCR423,
                              "Strongly agree" = 3, "Agree" = 2, "Slightly agree" = 1,
                              "Neither agree nor disagree" = 0, "Slightly disagree" = -1,
                              "Disagree" = -2, "Strongly disagree" = -3),
    
    policy_attitudes_index = rowMeans(dplyr::across(c(frei_att1, frei_att2, frei_att3)), na.rm = TRUE),
    
    # Treatment + state factor
    treat = ifelse(UCR419_text_split == "See UCR419 text", 1, 0),
    inputstate = as.factor(inputstate)
  )

# ============================================================
# 3. IPW weights to correct for ZIP missingness
# ============================================================

cces_full <- cces_data %>%
  mutate(has_zip = ifelse(!is.na(inputzip) & inputzip != "", 1, 0)) %>%
  filter(
    !is.na(imm_att_index),
    !is.na(econ_pers),
    !is.na(income_cont),
    !is.na(racial_resentment),
    !is.na(age),
    !is.na(gender_bin),
    !is.na(race_black),
    !is.na(race_asian),
    !is.na(race_hispanic),
    !is.na(race_other),
    !is.na(white),
    !is.na(educ_cont),
    !is.na(party_id),
    !is.na(homeowner),
    !is.na(urbancity),
    !is.na(region),
    !is.na(income_change),
    !is.na(price_change),
    !is.na(imm_inflation),
    !is.na(foreign_inflation),
    !is.na(employed)
  ) %>%
  mutate(rowid = row_number())

nrow(cces_full)  # check filtering

zip_missing_model <- glm(
  has_zip ~ gender_bin + age + white + educ_cont + income_cont +
    party_id + homeowner + employed + imm_att_index + income_change +
    price_change + imm_inflation + foreign_inflation + region + urbancity,
  data = cces_full,
  family = binomial(link = "logit")
)

summary(zip_missing_model)

# ---- GBM weights ----
weight_model_gbm <- WeightIt::weightit(
  has_zip ~ income_cont + racial_resentment + age + gender_bin +
    white + educ_cont + party_id + homeowner + urbancity + region +
    employed,
  data = cces_full,
  method = "gbm",
  estimand = "ATE"
)

cces_weighted <- cces_full %>%
  filter(has_zip == 1) %>%
  mutate(ipw_gbm_zip = weight_model_gbm$weights[cces_full$has_zip == 1])

summary(cces_weighted$ipw_gbm_zip)
mean(cces_weighted$ipw_gbm_zip > 5)
mean(cces_weighted$ipw_gbm_zip > 10)

# ---- Balance ----
bal_gbm <- cobalt::bal.tab(weight_model_gbm, un = TRUE)

core_vars <- c(
  "income_cont", "racial_resentment", "age", "gender_bin", "white",
  "educ_cont", "party_id", "homeowner", "urbancity", "region", "employed"
)

bal_gbm$Balance

love_plot_gbm <- cobalt::love.plot(
  bal_gbm,
  subset = core_vars,
  var.order = "unadjusted",
  abs = TRUE,
  threshold = 0.1,
  colors = c("grey60", "black"),
  stars = "raw"
)

print(love_plot_gbm)

# ============================================================
# 4. Merge Zillow + ACS data (Oct 2024 endpoints)
# ============================================================

# ---- Zillow (zillow1.csv): ZIP home values ----
zillow_wide <- readr::read_csv("zillow1.csv") %>%
  mutate(inputzip = stringr::str_pad(as.character(RegionName), 5, pad = "0"))

zillow_change <- zillow_wide %>%
  mutate(price_change_zip = 100 * (`2024-10-31` - `2019-10-31`) / `2019-10-31`) %>%
  dplyr::select(inputzip, price_change_zip)

zillow_latest <- zillow_wide %>%
  transmute(inputzip = inputzip, home_price_oct2024 = `2024-10-31`)

# ---- ACS median HH income (2019 vs 2024 ACS 5-year) ----
income_2019 <- tidycensus::get_acs(
  geography = "zip code tabulation area",
  variables = "B19013_001",
  year = 2019,
  survey = "acs5"
) %>%
  transmute(inputzip = str_pad(GEOID, 5, pad = "0"),
            income_2019 = estimate)

income_2024 <- tidycensus::get_acs(
  geography = "zip code tabulation area",
  variables = "B19013_001",
  year = 2024,
  survey = "acs5"
) %>%
  transmute(inputzip = str_pad(GEOID, 5, pad = "0"),
            income_2024 = estimate)

income_change <- income_2019 %>%
  left_join(income_2024, by = "inputzip") %>%
  mutate(income_pct_change_zip = 100 * (income_2024 - income_2019) / income_2019)

# ---- Merge into CCES (ZIP-observed sample) ----
merged_data <- cces_weighted %>%
  left_join(zillow_change, by = "inputzip") %>%
  left_join(income_change, by = "inputzip") %>%
  mutate(
    affordability_gap_zip = income_pct_change_zip - price_change_zip,
    price_to_income_ratio = price_change_zip / income_pct_change_zip
  ) %>%
  drop_na(affordability_gap_zip)

merged_data <- merged_data %>%
  left_join(zillow_latest, by = "inputzip") %>%
  left_join(income_2024 %>% rename(income_2024_static = income_2024), by = "inputzip") %>%
  mutate(static_affordability_ratio = home_price_oct2024 / income_2024_static) %>%
  drop_na(static_affordability_ratio)

# ============================================================
# 5. Standardize key continuous variables + blame variables
# ============================================================

merged_data_std <- merged_data %>%
  mutate(
    affordability_gap_zip_std      = scale(affordability_gap_zip)[, 1],
    price_to_income_ratio_std      = scale(price_to_income_ratio)[, 1],
    static_affordability_ratio_std = scale(static_affordability_ratio)[, 1],
    lat_ins_std                    = scale(lat_ins)[, 1],
    policy_attitudes_index_std     = scale(policy_attitudes_index)[, 1],
    educ_cont_std                  = scale(educ_cont)[, 1],
    income_cont_std                = scale(income_cont)[, 1],
    racial_resentment_std          = scale(racial_resentment)[, 1],
    imm_att_index_std              = scale(imm_att_index)[, 1],
    income_change_std              = scale(income_change)[, 1],
    
    blame_score = rowSums(across(c(imm_inflation, foreign_inflation)), na.rm = TRUE),
    blame_attribution = if_else(imm_inflation == 1 | foreign_inflation == 1, 1, 0, missing = NA_real_)
  )

# ============================================================
# 6. Summary statistics (weighted)
# ============================================================

vars_to_summarize <- c(
  "affordability_gap_zip", "affordability_gap_zip_std",
  "price_to_income_ratio", "price_to_income_ratio_std",
  "static_affordability_ratio", "static_affordability_ratio_std",
  "lat_ins", "lat_ins_std",
  "policy_attitudes_index", "policy_attitudes_index_std"
)

summary_data <- merged_data_std %>%
  dplyr::select(all_of(vars_to_summarize), ipw_gbm_zip) %>%
  drop_na()

weighted_stats <- function(x, w) {
  # NOTE: keep exactly your behavior, but ensure weights align to x non-missing
  ok <- !is.na(x) & !is.na(w)
  x <- x[ok]; w <- w[ok]
  c(
    min = min(x),
    median = as.numeric(Hmisc::wtd.quantile(x, weights = w, probs = 0.5)),
    mean = Hmisc::wtd.mean(x, weights = w),
    sd = sqrt(Hmisc::wtd.var(x, weights = w)),
    max = max(x)
  )
}

summary_table <- summary_data %>%
  dplyr::select(-ipw_gbm_zip) %>%
  purrr::map_dfr(~weighted_stats(.x, summary_data$ipw_gbm_zip), .id = "variable") %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))

summary_table

# ============================================================
# 7. Correlations
# ============================================================

key_vars <- merged_data_std %>%
  dplyr::select(
    ipw_gbm_zip,
    affordability_gap_zip_std,
    price_to_income_ratio_std,
    static_affordability_ratio_std,
    lat_ins_std,
    policy_attitudes_index_std
  )

cor_matrix <- cor(key_vars, use = "complete.obs")
print(round(cor_matrix, 3))

# Weight correlations with key vars (diagnostic)
summary(lm(lat_ins_std ~ ipw_gbm_zip, data = merged_data_std))
summary(lm(policy_attitudes_index_std ~ ipw_gbm_zip, data = merged_data_std))
summary(lm(treat ~ ipw_gbm_zip, data = merged_data_std))

summary(lm(affordability_gap_zip_std ~ ipw_gbm_zip, data = merged_data_std))
summary(lm(price_to_income_ratio_std ~ ipw_gbm_zip, data = merged_data_std))
summary(lm(static_affordability_ratio_std ~ ipw_gbm_zip, data = merged_data_std))

# Pairwise correlations among IVs + Holm correction
afford_vars <- merged_data_std %>%
  dplyr::select(
    affordability_gap_zip_std,
    price_to_income_ratio_std,
    static_affordability_ratio_std
  )

pairs <- combn(names(afford_vars), 2, simplify = FALSE)

tests <- lapply(pairs, function(p) {
  ct <- cor.test(afford_vars[[p[1]]], afford_vars[[p[2]]], use = "complete.obs")
  data.frame(
    pair = paste(p[1], "vs", p[2]),
    r = unname(ct$estimate),
    p = ct$p.value
  )
})

results <- dplyr::bind_rows(tests)
results$p_holm <- p.adjust(results$p, method = "holm")
results

# ============================================================
# 8. Direct effects models: lateral insecurity
# ============================================================

tmb_lat_unweighted_gap <- glmmTMB(
  lat_ins_std ~ affordability_gap_zip_std + treat + age_group + gender_bin + educ_cont_std +
    income_cont_std + employed + party_id + white + homeowner + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_lat_weighted_gap <- glmmTMB(
  lat_ins_std ~ affordability_gap_zip_std + treat + age_group + gender_bin + educ_cont_std +
    income_cont_std + employed + party_id + white + homeowner + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_lat_unweighted_ratio <- glmmTMB(
  lat_ins_std ~ price_to_income_ratio_std + treat + age_group + gender_bin + educ_cont_std +
    income_cont_std + employed + party_id + white + homeowner + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_lat_weighted_ratio <- glmmTMB(
  lat_ins_std ~ price_to_income_ratio_std + treat + age_group + gender_bin + educ_cont_std +
    income_cont_std + employed + party_id + white + homeowner + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~ price_to_income_ratio_std,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_lat_unweighted_static <- glmmTMB(
  lat_ins_std ~ static_affordability_ratio_std + treat + age_group + gender_bin + educ_cont_std +
    income_cont_std + employed + party_id + white + homeowner + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_lat_weighted_static <- glmmTMB(
  lat_ins_std ~ static_affordability_ratio_std + treat + age_group + gender_bin + educ_cont_std +
    income_cont_std + employed + party_id + white + homeowner + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

summary(tmb_lat_unweighted_gap); summary(tmb_lat_weighted_gap)
summary(tmb_lat_unweighted_ratio); summary(tmb_lat_weighted_ratio)
summary(tmb_lat_unweighted_static); summary(tmb_lat_weighted_static)

# ============================================================
# 9. Direct effects models: FREI policy preferences
# ============================================================

tmb_policy_unweighted_gap <- glmmTMB(
  policy_attitudes_index_std ~ affordability_gap_zip_std + treat + lat_ins_std + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + homeowner + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_policy_weighted_gap <- glmmTMB(
  policy_attitudes_index_std ~ affordability_gap_zip_std + treat + lat_ins_std + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + homeowner + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_policy_unweighted_ratio <- glmmTMB(
  policy_attitudes_index_std ~ price_to_income_ratio_std + treat + lat_ins_std + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + homeowner + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_policy_weighted_ratio <- glmmTMB(
  policy_attitudes_index_std ~ price_to_income_ratio_std + treat + lat_ins_std + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + homeowner + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_policy_unweighted_static <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std + treat + lat_ins_std + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + homeowner + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_policy_weighted_static <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std + treat + lat_ins_std + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + homeowner + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

summary(tmb_policy_unweighted_gap); summary(tmb_policy_weighted_gap)
summary(tmb_policy_unweighted_ratio); summary(tmb_policy_weighted_ratio)
summary(tmb_policy_unweighted_static); summary(tmb_policy_weighted_static)

# ============================================================
# 10. Moderation by housing tenure (homeowner): lat insecurity
# ============================================================

tmb_lat_interact_unweighted_gap <- glmmTMB(
  lat_ins_std ~ affordability_gap_zip_std * homeowner + treat + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_lat_interact_weighted_gap <- glmmTMB(
  lat_ins_std ~ affordability_gap_zip_std * homeowner + treat + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_lat_interact_unweighted_ratio <- glmmTMB(
  lat_ins_std ~ price_to_income_ratio_std * homeowner + treat + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_lat_interact_weighted_ratio <- glmmTMB(
  lat_ins_std ~ price_to_income_ratio_std * homeowner + treat + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_lat_interact_unweighted_static <- glmmTMB(
  lat_ins_std ~ static_affordability_ratio_std * homeowner + treat + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_lat_interact_weighted_static <- glmmTMB(
  lat_ins_std ~ static_affordability_ratio_std * homeowner + treat + age_group +
    gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

summary(tmb_lat_interact_unweighted_gap); summary(tmb_lat_interact_weighted_gap)
summary(tmb_lat_interact_unweighted_ratio); summary(tmb_lat_interact_weighted_ratio)
summary(tmb_lat_interact_unweighted_static); summary(tmb_lat_interact_weighted_static)

# ============================================================
# 11. Moderation by housing tenure (homeowner): FREI policy
# ============================================================

tmb_frei_interact_unweighted_gap <- glmmTMB(
  policy_attitudes_index_std ~ affordability_gap_zip_std * homeowner + treat + lat_ins_std +
    age_group + gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_frei_interact_weighted_gap <- glmmTMB(
  policy_attitudes_index_std ~ affordability_gap_zip_std * homeowner + treat + lat_ins_std +
    age_group + gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_frei_interact_unweighted_ratio <- glmmTMB(
  policy_attitudes_index_std ~ price_to_income_ratio_std * homeowner + treat + lat_ins_std +
    age_group + gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_frei_interact_weighted_ratio <- glmmTMB(
  policy_attitudes_index_std ~ price_to_income_ratio_std * homeowner + treat + lat_ins_std +
    age_group + gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

tmb_frei_interact_unweighted_static <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std * homeowner + treat + lat_ins_std +
    age_group + gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_frei_interact_weighted_static <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std * homeowner + treat + lat_ins_std +
    age_group + gender_bin + educ_cont_std + income_cont_std + employed +
    party_id + white + racial_resentment_std + income_change_std +
    urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

summary(tmb_frei_interact_unweighted_gap); summary(tmb_frei_interact_weighted_gap)
summary(tmb_frei_interact_unweighted_ratio); summary(tmb_frei_interact_weighted_ratio)
summary(tmb_frei_interact_unweighted_static); summary(tmb_frei_interact_weighted_static)

# ============================================================
# 12. Appendix C: Static affordability supplemental tests
# ============================================================

merged_data_std <- merged_data_std %>%
  mutate(
    recent_unaffordability = ifelse(affordability_gap_zip < -5, 1, 0),
    recent_unaffordability = factor(recent_unaffordability, levels = c(0, 1),
                                    labels = c("No", "Yes")),
    homeowner = as.numeric(homeowner)
  )

# ---- Font handling ----
# NOTE: this is macOS-specific.
font_add("Times New Roman", regular = "/System/Library/Fonts/Supplemental/Times New Roman.ttf")
showtext_auto()

# ---- Distribution figure ----
aff_gap_plot <- ggplot(merged_data_std, aes(x = affordability_gap_zip)) +
  geom_histogram(bins = 30, fill = "grey80", color = "black", linewidth = 0.3) +
  geom_vline(xintercept = -5, color = "red", linetype = "dashed", linewidth = 0.8) +
  labs(x = "Affordability Gap", y = "Number of Respondents") +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    axis.title.x = element_text(size = 11, face = "bold"),
    axis.title.y = element_text(size = 11, face = "bold"),
    axis.text = element_text(size = 10, color = "black"),
    panel.grid.major = element_line(linewidth = 0.2),
    panel.grid.minor = element_blank()
  )

ggsave("aff_gap_histogram.pdf", aff_gap_plot,
       width = 6, height = 4, dpi = 300, device = cairo_pdf)

# ---- Model 1: Static × Recent unaffordability ----
tmb_policy_recent_vs_static_un <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std * recent_unaffordability +
    treat + lat_ins_std + age_group + gender_bin + educ_cont_std + income_cont_std +
    employed + party_id + white + homeowner + racial_resentment_std +
    income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_policy_recent_vs_static <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std * recent_unaffordability +
    treat + lat_ins_std + age_group + gender_bin + educ_cont_std + income_cont_std +
    employed + party_id + white + homeowner + racial_resentment_std +
    income_change_std + urbancity + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

summary(tmb_policy_recent_vs_static_un)
summary(tmb_policy_recent_vs_static)

mfx <- ggeffects::ggpredict(
  tmb_policy_recent_vs_static,
  terms = c("static_affordability_ratio_std", "recent_unaffordability")
) %>% as.data.frame()

policy_plot <- ggplot(mfx, aes(x = x, y = predicted, color = group)) +
  geom_line(linewidth = 0.9) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group),
              alpha = 0.25, color = NA) +
  labs(
    x = "Static Affordability Ratio (Standardized)",
    y = "Predicted FREI Policy Attitudes",
    color = "Recent Unaffordability",
    fill = "Recent Unaffordability"
  ) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 11)
  )

ggsave("ME_plot1.pdf", policy_plot,
       width = 6.5, height = 4.5, units = "in", dpi = 300, device = cairo_pdf)

# ---- Model 2: Static × Urbanicity ----
tmb_policy_urbancity_vs_static_un <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std * urbancity +
    treat + lat_ins_std + age_group + gender_bin + educ_cont_std + income_cont_std +
    employed + party_id + white + homeowner + racial_resentment_std +
    income_change_std + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std
)

tmb_policy_urbancity_vs_static <- glmmTMB(
  policy_attitudes_index_std ~ static_affordability_ratio_std * urbancity +
    treat + lat_ins_std + age_group + gender_bin + educ_cont_std + income_cont_std +
    employed + party_id + white + homeowner + racial_resentment_std +
    income_change_std + (1 | inputstate),
  dispformula = ~1,
  data = merged_data_std,
  weights = ipw_gbm_zip
)

summary(tmb_policy_urbancity_vs_static_un)
summary(tmb_policy_urbancity_vs_static)

mfx_urb <- ggeffects::ggpredict(
  tmb_policy_urbancity_vs_static,
  terms = c("static_affordability_ratio_std", "urbancity")
) %>% as.data.frame()

urb_plot_faceted <- ggplot(mfx_urb, aes(x = x, y = predicted)) +
  geom_line(aes(color = group), linewidth = 0.9) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group),
              alpha = 0.25, color = NA) +
  facet_wrap(~ group, ncol = 2, scales = "free_y") +
  labs(
    x = "Static Affordability Ratio (Standardized)",
    y = "Predicted FREI Policy Attitudes"
  ) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    strip.text = element_text(size = 11, face = "bold"),
    axis.title.x = element_text(size = 12, margin = margin(t = 10)),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text = element_text(size = 10),
    panel.spacing = grid::unit(1, "lines"),
    legend.position = "none"
  )

ggsave("ME_plot2.pdf", urb_plot_faceted,
       width = 7, height = 5.5, units = "in", dpi = 300, device = cairo_pdf)

# ============================================================
# 13. Main effects plots + interaction plots + tenure moderation plots
# ============================================================

# -----------------------------
# Plot interaction coefficients (Affordability × Treatment)
# -----------------------------

library(ggplot2)
library(tibble)
library(dplyr)
library(showtext)

# Set font for macOS
font_add("Times New Roman", regular = "/System/Library/Fonts/Supplemental/Times New Roman.ttf")
showtext_auto()

# Data for interaction coefficients
interaction_data <- tribble(
  ~model, ~type, ~estimate, ~se,
  "Affordability Gap", "Unweighted",  0.08084, 0.13562,
  "Affordability Gap", "Weighted",    0.10397, 0.08430,
  "Price-to-Income",   "Unweighted", -0.07896, 0.29359,
  "Price-to-Income",   "Weighted",   -0.06257, 0.17407,
  "Static Ratio",      "Unweighted",  0.38403, 0.24160,
  "Static Ratio",      "Weighted",    0.40124, 0.15382
) %>%
  mutate(
    model = factor(model, levels = c("Affordability Gap", "Price-to-Income", "Static Ratio")),
    model_type = paste0(model, " (", type, ")"),
    model_type = factor(model_type, levels = rev(unique(model_type[order(model)]))),
    lower = estimate - 1.96 * se,
    upper = estimate + 1.96 * se
  )

# Save PDF
pdf("interaction_plot.pdf", width = 6, height = 4, family = "Times New Roman")

ggplot(interaction_data, aes(x = estimate, y = model_type, color = model)) +
  geom_point(size = 2) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray40") +
  xlab("Interaction Term Coefficient (Affordability × Treatment)") +
  ylab(NULL) +
  theme_minimal(base_size = 10, base_family = "Times New Roman") +
  theme(
    axis.text.y = element_text(size = 10),
    plot.title = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 10))
  )

dev.off()

# -----------------------------
# Marginal effects plot (Static affordability × Treatment)
# -----------------------------

library(ggplot2)
library(ggeffects)
library(showtext)

# Load font
font_add("Times New Roman", regular = "/System/Library/Fonts/Supplemental/Times New Roman.ttf")
showtext_auto()

# Get marginal effects for interaction
me_static <- ggpredict(int_policy_weighted_static, terms = c("static_affordability_ratio_std", "treat"))

# Create marginal effects plot with pastel red and blue
pdf("ME_plot_static.pdf", width = 6, height = 4, family = "Times New Roman")

ggplot(me_static, aes(x = x, y = predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.3, color = NA) +
  geom_line(size = 1) +
  scale_color_manual(values = c("0" = "#92c5de", "1" = "#f4a582")) +  # blue, red
  scale_fill_manual(values = c("0" = "#92c5de", "1" = "#f4a582")) +  # blue, red
  labs(
    x = "Static Affordability Ratio (standardized)",
    y = "Predicted Restrictive FREI Preferences (standardized)",
    color = "Treatment",
    fill = "Treatment"
  ) +
  theme_minimal(base_size = 10, base_family = "Times New Roman") +
  theme(
    legend.position = "top",
    axis.title.y = element_text(margin = margin(r = 10)),
    axis.title.x = element_text(margin = margin(t = 10))
  )

dev.off()

# -----------------------------
# Plot moderation by housing tenure (Homeowner vs Non-Homeowner)
# -----------------------------

library(ggeffects)
library(ggplot2)

# LATERAL INSECURITY PREDICTIONS
gap_latins <- ggpredict(tmb_lat_interact_weighted_gap, terms = c("affordability_gap_zip_std", "homeowner"))
ratio_latins <- ggpredict(tmb_lat_interact_weighted_ratio, terms = c("price_to_income_ratio_std", "homeowner"))
static_latins <- ggpredict(tmb_lat_interact_weighted_static, terms = c("static_affordability_ratio_std", "homeowner"))

# FREI POLICY ATTITUDE PREDICTIONS
gap_frei <- ggpredict(tmb_frei_interact_weighted_gap, terms = c("affordability_gap_zip_std", "homeowner"))
ratio_frei <- ggpredict(tmb_frei_interact_weighted_ratio, terms = c("price_to_income_ratio_std", "homeowner"))
static_frei <- ggpredict(tmb_frei_interact_weighted_static, terms = c("static_affordability_ratio_std", "homeowner"))

plot_interaction <- function(pred_data, title, ylab) {
  ggplot(pred_data, aes(x = x, y = predicted, color = group, fill = group)) +
    geom_line(linewidth = 1) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, linetype = 0) +
    labs(title = title, x = "Affordability (Standardized)", y = ylab) +
    theme_minimal(base_family = "Times New Roman", base_size = 13) +
    scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("Non-Homeowner", "Homeowner")) +
    scale_fill_manual(values = c("#0072B2", "#D55E00"), labels = c("Non-Homeowner", "Homeowner")) +
    theme(
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.text = element_text(size = 13),
      axis.title = element_text(size = 14),
      axis.text = element_text(size = 12),
      plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
    )
}

# Row 1: Lateral Insecurity
p_gap_latins    <- plot_interaction(gap_latins,    "(a) Affordability Gap", "Lateral Insecurity")
p_ratio_latins  <- plot_interaction(ratio_latins,  "(b) Price-to-Income Ratio", "Lateral Insecurity")
p_static_latins <- plot_interaction(static_latins, "(c) Static Affordability", "Lateral Insecurity")

# Row 2: FREI Attitudes
p_gap_frei    <- plot_interaction(gap_frei,    "(d) Affordability Gap", "Restrictive FREI Attitudes")
p_ratio_frei  <- plot_interaction(ratio_frei,  "(e) Price-to-Income Ratio", "Restrictive FREI Attitudes")
p_static_frei <- plot_interaction(static_frei, "(f) Static Affordability", "Restrictive FREI Attitudes")

library(patchwork)

final_plot <- (p_gap_latins | p_ratio_latins | p_static_latins) /
  (p_gap_frei | p_ratio_frei | p_static_frei) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

# Display
print(final_plot)

library(extrafont)

# Only do this once:
# font_import(pattern = "Times", prompt = FALSE)
# loadfonts(device = "pdf")  # Run this once per R session

# Save plot with Times New Roman
ggsave("interaction_plots_combined.pdf", final_plot, width = 11, height = 7.5, family = "Times New Roman")

# -----------------------------
# Plot main (non-interaction) models: weighted coefficients
# -----------------------------

library(glmmTMB)
library(broom.mixed)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(extrafont)

# If running for the first time:
# font_import()  # Only once
loadfonts(device = "pdf")

# Define consistent color scheme
afford_colors <- c(
  "Afford. Gap" = "#1f77b4",
  "Price-to-Income Ratio" = "#2ca02c",
  "Static Ratio" = "#d62728"
)

# Extractor function
get_main_effect <- function(model, term_label, model_label) {
  broom.mixed::tidy(model, effects = "fixed") %>%
    filter(term == term_label) %>%
    mutate(model = model_label)
}

# Extract data
lat_data <- bind_rows(
  get_main_effect(tmb_lat_weighted_gap, "affordability_gap_zip_std", "Afford. Gap"),
  get_main_effect(tmb_lat_weighted_ratio, "price_to_income_ratio_std", "Price-to-Income Ratio"),
  get_main_effect(tmb_lat_weighted_static, "static_affordability_ratio_std", "Static Ratio")
)

policy_data <- bind_rows(
  get_main_effect(tmb_policy_weighted_gap, "affordability_gap_zip_std", "Afford. Gap"),
  get_main_effect(tmb_policy_weighted_ratio, "price_to_income_ratio_std", "Price-to-Income Ratio"),
  get_main_effect(tmb_policy_weighted_static, "static_affordability_ratio_std", "Static Ratio")
)

# Reorder IVs for consistent top-to-bottom order
order_levels <- c("Static Ratio", "Price-to-Income Ratio", "Afford. Gap")

lat_data$model <- factor(lat_data$model, levels = rev(order_levels))
policy_data$model <- factor(policy_data$model, levels = rev(order_levels))

# Set base theme with larger y-axis labels
theme_set(
  theme_minimal(base_size = 12, base_family = "Times New Roman") +
    theme(
      axis.text.y = element_text(size = 13),
      plot.title = element_text(face = "bold", size = 14)
    )
)

# Lateral Insecurity Plot
plot_lat <- ggplot(lat_data, aes(x = estimate, y = model, color = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray30", linewidth = 0.6) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error,
                     xmax = estimate + 1.96 * std.error),
                 height = 0.2, linewidth = 0.8) +
  scale_color_manual(values = afford_colors) +
  labs(
    title = "Lateral Insecurity",
    x = "Coefficient Estimate", y = NULL, color = NULL
  )

# Policy Attitudes Plot
plot_policy <- ggplot(policy_data, aes(x = estimate, y = model, color = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray30", linewidth = 0.6) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error,
                     xmax = estimate + 1.96 * std.error),
                 height = 0.2, linewidth = 0.8) +
  scale_color_manual(values = afford_colors) +
  labs(
    title = "FREI Policy Preferences",
    x = "Coefficient Estimate", y = NULL, color = NULL
  )

# Combine vertically
combined_plot <- ggarrange(plot_lat, plot_policy,
                           ncol = 1, nrow = 2,
                           common.legend = TRUE,
                           legend = "bottom")

# Save to PDF
ggsave("main_effects_plot.pdf", combined_plot,
       width = 7.5, height = 10, dpi = 300)


# ============================================================
# 14. Dispersion tests (LRT) + random intercept tests
# ============================================================

# Dispersion tests for lateral insecurity (weighted models)
m0 <- update(tmb_lat_weighted_gap, dispformula = ~1)
m1 <- update(tmb_lat_weighted_gap, dispformula = ~ affordability_gap_zip_std)
anova(m0, m1)

m0 <- update(tmb_lat_weighted_ratio, dispformula = ~1)
m1 <- update(tmb_lat_weighted_ratio, dispformula = ~ price_to_income_ratio_std)
anova(m0, m1)

m0 <- update(tmb_lat_weighted_static, dispformula = ~1)
m1 <- update(tmb_lat_weighted_static, dispformula = ~ static_affordability_ratio_std)
anova(m0, m1)

# Dispersion tests for policy attitudes (weighted models)
m0 <- update(tmb_policy_weighted_gap, dispformula = ~1)
m1 <- update(tmb_policy_weighted_gap, dispformula = ~ affordability_gap_zip_std)
anova(m0, m1)

m0 <- update(tmb_policy_weighted_ratio, dispformula = ~1)
m1 <- update(tmb_policy_weighted_ratio, dispformula = ~ price_to_income_ratio_std)
anova(m0, m1)

m0 <- update(tmb_policy_weighted_static, dispformula = ~1)
m1 <- update(tmb_policy_weighted_static, dispformula = ~ static_affordability_ratio_std)
anova(m0, m1)

# Random intercept test helper
lr_random_intercept_test <- function(model_re, data, weights, label) {
  model_fe <- update(model_re, . ~ . - (1 | inputstate))
  model_fe_refit <- update(model_fe, data = data, weights = weights)
  model_re_refit <- update(model_re, data = data, weights = weights)
  lr <- anova(model_fe_refit, model_re_refit)
  
  tibble::tibble(
    model = label,
    df_fe = lr$Df[1],
    df_re = lr$Df[2],
    chisq = lr$Chisq[2],
    df = lr$`Chi Df`[2],
    p = lr$`Pr(>Chisq)`[2],
    aic_fe = lr$AIC[1],
    aic_re = lr$AIC[2]
  )
}

lr_table <- dplyr::bind_rows(
  lr_random_intercept_test(tmb_lat_weighted_gap,    merged_data_std, merged_data_std$ipw_gbm_zip, "Lat. Ins. -- Aff. Gap"),
  lr_random_intercept_test(tmb_lat_weighted_ratio,  merged_data_std, merged_data_std$ipw_gbm_zip, "Lat. Ins. -- Ratio"),
  lr_random_intercept_test(tmb_lat_weighted_static, merged_data_std, merged_data_std$ipw_gbm_zip, "Lat. Ins. -- Static"),
  lr_random_intercept_test(tmb_policy_weighted_gap,    merged_data_std, merged_data_std$ipw_gbm_zip, "Policy Pref. -- Aff. Gap"),
  lr_random_intercept_test(tmb_policy_weighted_ratio,  merged_data_std, merged_data_std$ipw_gbm_zip, "Policy Pref. -- Ratio"),
  lr_random_intercept_test(tmb_policy_weighted_static, merged_data_std, merged_data_std$ipw_gbm_zip, "Policy Pref. -- Static")
)

lr_table
